

%Run programs
clc
clear all
close all
% 
%Initial steady-state
run Parameters
run BGP_growth

ss_c=C;
ss_y=Y;
ss_yT=YT;
ss_yNT=YNT;
ss_ha=Ha;
ss_hr=H_r;
ss_g=g(length(g));
ss_TT=Tnew;
ss_TNT=TNT;
ss_pi=Pi;
ss_piT=PiT;
ss_piNT=PiNT;
ss_p=P;
ss_pT=PT;
ss_pNT=PNT;
ss_PishareT=pi_shareT;
ss_PishareNT=pi_shareNT;
ss_w=w;
ss_vinnovator=vinnovator;
ss_jinnovator=Jinnovator;
ss_v=vadopter;
ss_j=Jadopter;
ss_vinnov=vinnov;
ss_chi = chi;
ss_YL=Y_L;
ss_eps=eps;
ss_Z=ZZ;
ss_AT=A_Z;
ss_rp=roy;
ss_IM=LHS;
ss_EX=RHS;
ss_IMT=IM_T;
ss_EXT=EX_T;
ss_IMNT=IM_NT;
ss_EXNT=EX_NT;
ss_Yw=Yw;
ss_Deltaeps=Delta_eps;
ss_rhochi = rho_chi; 

save Initial_values.mat ss_rhochi ss_chi ss_Yw Q ss_TNT ss_y ss_yT ss_yNT ss_pT ss_pNT ss_IMT ss_EXT ss_IMNT ss_EXNT ss_Deltaeps ss_c ss_y ss_ha ss_hr ss_g ss_TT ss_pi ss_piT ss_piNT ss_p ss_PishareT ss_PishareNT ss_w ss_vinnovator ss_jinnovator ss_v ss_j ss_vinnov ss_YL ss_eps ss_Z ss_AT ss_rp ss_EX ss_IM;
save Q.mat Q
save epsbar.mat epsbar
save lam.mat lam


%%Counterfactual steady-state
run Parameters_counterf
run BGP_growth_counterf


ss1_c=C;
ss1_y=Y;
ss1_yT=YT;
ss1_yNT=YNT;
ss1_ha=Ha;
ss1_hr=H_r;
ss1_g=g(length(g));
ss1_TT=Tnew;
ss1_TNT=TNT;
ss1_pi=Pi;
ss1_piT=PiT;
ss1_piNT=PiNT;
ss1_p=P;
ss1_pT=PT;
ss1_pNT=PNT;
ss1_PishareT=pi_shareT;
ss1_PishareNT=pi_shareNT;
ss1_w=w;
ss1_vinnovator=vinnovator;
ss1_chi = chi;
ss1_jinnovator=Jinnovator;
ss1_v=vadopter;
ss1_j=Jadopter;
ss1_vinnov=vinnov;
ss1_YL=Y_L;
ss1_eps=eps;
ss1_Z=ZZ;
ss1_AT=A_Z;
ss1_rp=roy;
ss1_IM=LHS;
ss1_EX=RHS;
ss1_IMT=IM_T;
ss1_EXT=EX_T;
ss1_IMNT=IM_NT;
ss1_EXNT=EX_NT;
ss1_Yw=Yw;
ss1_Deltaeps=Delta_eps;
ss1_rhochi = rho_chi; 

save Final_values.mat ss1_rhochi ss1_chi ss1_Yw ss1_TNT ss1_y ss1_yT ss1_yNT ss1_pT ss1_pNT ss1_IMT ss1_EXT ss1_IMNT ss1_EXNT ss1_Deltaeps ss1_c ss1_y ss1_ha ss1_hr ss1_g ss1_TT ss1_pi ss1_piT ss1_piNT ss1_p ss1_PishareT ss1_PishareNT ss1_w ss1_vinnovator ss1_jinnovator ss1_v ss1_j ss1_vinnov ss1_YL ss1_eps ss1_Z ss1_AT ss1_rp ss1_EX ss1_IM;
save epsbar.mat epsbar
% Nash_equilibrium.m

% Clear all variables and close all figures
clear all;

% Number of iterations
num_iterations = 100;

% Initial guesses for the equilibrium strategy parameters
initial_guess_us = 0; 
initial_guess_chn = [1.2, 1.2]; 

% Initialize parameters for the first iteration
optimal_params_us = initial_guess_us;
optimal_params_chn = initial_guess_chn;

% Initialize previous parameters for convergence check
previous_params_us = optimal_params_us;
previous_params_chn = optimal_params_chn;

% Loop for iterations
for iter = 1:num_iterations
    % Run Dynare script to simulate the model and save variables for the US
    dynare 'Dynamic_GFT_exponents.mod';
    save DynamicGFT_US.mat resultsgc1 resultsC1 resultsgc3 resultsC3;

   
    % Load necessary variables from the MAT files for the US
    load DynamicGFT_US.mat;
    
    % Define the objective function for the Nash equilibrium for the US
    objective_us = @(params_us) -compute_welfare(params_us, optimal_params_chn, resultsgc1, resultsC1, resultsgc3, resultsC3);

    % Solve for the Nash equilibrium for the US
    [optimal_params_us, max_welfare_us] = fminsearch(objective_us, initial_guess_us);
 % Run Dynare script to simulate the model and save variables for China
    dynare 'Dynamic_GFT_exponents.mod';
    save DynamicGFT_CHN.mat resultsgc1 resultsC1 resultsgc3 resultsC3;

    % Load necessary variables from the MAT files for China
    load DynamicGFT_CHN.mat;
    
    % Define the objective function for the Nash equilibrium for China
    objective_chn = @(params_chn) -compute_welfare(optimal_params_us, params_chn, resultsgc1, resultsC1, resultsgc3, resultsC3);

    % Solve for the Nash equilibrium for China
    [optimal_params_chn, max_welfare_chn] = fminsearch(objective_chn, initial_guess_chn);

    % Display results for the current iteration
    disp(['Iteration ', num2str(iter)]);
    disp('Optimal Parameters for the US:');
    disp(optimal_params_us);
    disp(['Max Welfare for the US: ', num2str(-max_welfare_us)]);
    disp('Optimal Parameters for China:');
    disp(optimal_params_chn);
    disp(['Max Welfare for China: ', num2str(-max_welfare_chn)]);

    % Check for convergence
    % For example, you could check if the changes in parameters are small enough
    convergence_threshold = 1e-6;
    if max(abs(optimal_params_us - previous_params_us)) < convergence_threshold && ...
        max(abs(optimal_params_chn - previous_params_chn)) < convergence_threshold
        disp('Converged to a fixed point.');
        break;
    end

    % Update previous parameters for the next iteration
    previous_params_us = optimal_params_us;
    previous_params_chn = optimal_params_chn;
end

% Function to compute welfare based on policy functions
function welfare = compute_welfare(params_us, params_chn, resultsgc1, resultsC1, resultsgc3, resultsC3)
    tau_US = params_us(1);
    
    IP_CHNUSA = params_chn(1);
    IP_CHNCHN = params_chn(2);

    % Use resultsgc1, resultsC1, resultsgc3, and resultsC3 as needed
    TT = 150;
    bet = 0.985;
    
    Uold1 = sum((bet.^([1:TT]).*log((1+resultsgc1(1)).^[1:TT])));
    prodgc1 = cumprod((1+resultsgc1(1:TT)));
    Unew1 = sum((bet.^([1:TT]).*log(prodgc1)'));
    GFT1 = [exp((1-bet)*((log(exp(resultsgc1(TT))*(1+resultsgc1(TT)))-log(exp(resultsC1(1))*(1+resultsgc1(1))))/(1-bet)+Unew1-Uold1))-1]*100;

    % Calculate and return welfare based on model outcomes for the US
    welfare_us = sum(GFT1);  % Sum over the vector to obtain a scalar
    
    Uold3 = sum((bet.^([1:TT]).*log((1+resultsgc3(1)).^[1:TT])));
    prodgc3 = cumprod((1+resultsgc3(1:TT)));
    Unew3 = sum((bet.^([1:TT]).*log(prodgc3)'));
    GFT3 = [exp((1-bet)*((log(exp(resultsgc3(TT))*(1+resultsgc3(TT)))-log(exp(resultsC3(1))*(1+resultsgc3(1))))/(1-bet)+Unew3-Uold3))-1]*100;

    % Calculate and return welfare based on model outcomes for China
    welfare_chn = sum(GFT3);  % Sum over the vector to obtain a scalar

    % Combine individual welfare functions (you might use a different aggregation method)
    welfare = welfare_us + welfare_chn;
end
